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We present three-dimensional numerical simulations of turbulent combustion converting a neutron 
star into a quark star. Hadronic matter, described by a micro-physical finite-temperature equation 
of state, is converted into strange quark matter. We assume this phase, represented by a bag-model 
equation of state, to be absolutely stable. Following the example of thermonuclear burning in white 
dwarfs leading to Type la supernovae, we treat the conversion process as a potentially turbulent 
deflagration. Solving the non-relativistic Euler equations using established numerical methods we 
conduct large eddy simulations including an elaborate subgrid scale model, while the propagation 
of the conversion front is modeled with a level-set method. Our results show that for large parts 
of the parameter space the conversion becomes turbulent and therefore significantly faster than in 
the laminar case. Despite assuming absolutely stable strange quark matter, in our hydrodynamic 
approximation an outer layer remains in the hadronic phase, because the conversion front stops 
when it reaches conditions under which the combustion is no longer exothermic. 



I. INTRODUCTION 

Based on earlier work by Bodmer [l[ and Itoh Wit- 
ten Q suggested that a mixture of about the same num- 
ber of u-, d- and s-quarks, called strange quark matter 
(SQM), was the true ground state of matter, whereas or- 
dinary nuclear matter is only a metastable, yet usually 
extremely long-lived state. This conjecture, known today 
as strange matter hypothesis, was discussed lively ever 
since, but no final verdict about its correctness could be 
made because the equation of state (EoS) of cold dense 
matter is still largely unknown. Matter in this extreme 
state is inaccessible to laboratory experiments; compact 
stars, however, offer a possibility to test the strange mat- 
ter hypothesis. Shortly after Witten's work also Haensel 
et al. [ ]1 and Alcock et al. Q proposed strange stars, 
compact stars consisting entirely of SQM. Alcock et al. 
Q based their work on the idea that compact stars are 
not born as strange stars, but as hadronic neutron stars, 
which later are converted into strange stars or hybrid 
stars - compact stars consisting of a quark core and 
hadronic outer layers. 

Hadronic matter does not decay into SQM sponta- 
neously, even though it would be energetically favorable, 
because this process would require a large amount of si- 
multaneous weak reactions - the probability for this to 
happen is vanishingly low. But if some SQM already ex- 
ists inside a neutron star, the diffusion of s-quarks from 
this seed into the surrounding hadronic matter would 
convert it into SQM. This conversion process should take 
place in a confined region and on length-scales small com- 
pared to the size of the star. It is expected to occur 
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only if the conversion releases energy, that is, if it is an 
exothermic process. The described situation is therefore 
similar to the propagation of a chemical flame, or even 
more similar to the thermonuclear burning inside a white 
dwarf during a Type la supernova (SN la). Thus, it is 
natural to think of the conversion of hadronic matter into 
SQM as a "combustion" . In the spirit of this analogy, we 
will sometimes refer to the conversion as "burning" and 
to the conversion front as "flame front" . Alcock et al. [j| 
were the first to suggest that a strange star may originate 
from a combustion of an ordinary neutron star. They also 
considered how a SQM seed which subsequently triggers 
the conversion into a strange star may come about and 
described various possibilities by either internal nucle- 
ation or external seeding. Subsequently the idea of a 
combustion was discussed in more detail by various au- 
thors @-[3. 

The laminar conversion velocity was first estimated 
by Olinto [l4|, and, with similar results, by Heiselberg 
et al. [lj|. Based on their results, Olesen and Madsen 
calculated the burning of a neutron star using a one- 
dimensional model with laminar burning and obtained 
conversion timescales from 10 _1 s to 10 2 s. Horvath and 
Benvenuto [(| suggested that the combustion should be 
turbulent due to various instabilities of the conversion 
front and therefore the conversion velocity should be en- 
hanced considerably (see [l6| for a recent update). Lu- 
gones et al. [H and Lugones and Benvenuto |ldl | pointed 
out the importance of the conditions for an exothermic 
combustion. The combustion mode was discussed from 
a hydrodynamic point of view also by Cho et al. 
Tokareva and Nusser [ll[ and Drago et al. where 
the latter expected the burning to be subsonic, although 
accelerated by turbulence. New ideas concerning the ini- 
tial seeding were recently published by Perez-Garcia et al. 
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[17j . They suggested that the self-annihilation of weakly 
interacting dark matter particles (WIMPs) inside a neu- 
tron star may provide a SQM seed. Recently, hydro- 
dynamic simulations of the combustion front were pre- 
sented by Niebergal et al. [IH . Their results of the lami- 
nar conversion velocity differed strongly from earlier esti- 
mates. On the observational side Leahy and Ouyed [l8j . 
extended in Ouyed et al. fl9| . examined the supernova 
SN 2006gy and suggested that this extremely luminous 
event can be explained by a "quark nova" - the transition 
of the newly formed neutron star to a strange quark star 
shortly after a core collapse supernova of a very massive 
star. 

Here we study the dynamical behavior of the conver- 
sion inside a neutron star. We model the conversion as a 
combustion, particularly as a subsonic deflagration. As 
mentioned above, it is widely assumed that the conver- 
sion process turns turbulent [e.g. 0, [IH, [l6[, but dynam- 
ical, multi-dimensional simulations have never been per- 
formed. Thus, our main focus will be to explore if and 
how turbulent motion occurs during the conversion pro- 
cess and to which consequences for the final state of the 
neutron star this may lead. 

This work is organized as follows: In Section Ql] we de- 
scribe the EoS that we use in our calculations. In Section 
IIIII we introduce our concept of modeling the conversion 
as a turbulent combustion, and in Section IIVI our nu- 
merical method is explained. We present numerical sim- 
ulations and their results in Section |V] and conclude in 
Section ED 



II. EQUATION OF STATE 

A. Equation of State for Hadronic Matter 

We consider the two micro-physical, finite temperature 
EoS which are most frequently used in simulations of 
astrophysical events such as core collapse supernovae and 
neutron star mergers: the EoS by Lattimer and Swesty 
20] (LS EoS) and by Shen et al. [H[ (Shen EoS). The 
LS EoS is based on a non-relativistic liquid drop model 
with an incompressibility modulus of K — 180 MeV. For 
calculating the Shen EoS relativistic mean field theory 
was applied, here K = 280 MeV is adopted. 

The recent measurement of the Shapiro delay of the 
binary millisecond pulsar J1614-2230 [22j yields a grav- 
itational mass of the pulsar of M = (1.97 ± 0.04) M Q . 
In contrast to the Shen EoS, the LS EoS is rather soft. 
Consequently it leads to a maximum mass for a hadronic 
non-rotating neutron star of only M^| x ~ 1.8-Mq and 
is therefore in conflict with the observation of pulsar 
J1614-2230. We nevertheless use the LS EoS in this work, 
because we do not claim to conduct realistic simulations 
but we rather see our work as a first step into this so 
far mostly unexplored field. An alternative would be to 
change the incompressibility modulus of the LS EoS to 
K = 220 MeV, which leads to a maximum mass com- 



patible with the observations. We discuss this possibility 
briefly in Section IHI Bl 

For simplicity we assume for all our calculations a con- 
stant low proton fraction Y p . Variations of its value, par- 
ticularly assuming /3-equilibrium, do not lead to a sig- 
nificant change of our results, as is shown exemplary in 
Section IIIIBI In the same section we explain that for 
physical reasons it turned out that it was impossible to 
use the Shen EoS, thus we perform all our simulations 
using the LS EoS. 



B. Equation of State for Strange Quark Matter 

We describe SQM by a simple bag model for finite tem- 
(23 



peratures [23JJ based on the MIT bag model (2J. This 
model treats SQM as ideal Fermi gases of massless and 
non-interacting u-, d-, and s-quarks inside a confining 
bag. Since in this approximation the quarks can be de- 
scribed by only one chemical potential, the resulting an- 
alytic expressions for pressure P, energy density e and 
baryon number density n as function of the bag constant 
B and the chemical potential /x are [ill, HE] 



P 



19 

36 
19 



TT 2 T 4 



12 

tV 



TT 2 T 4 



1 



TV 2 
T 2 /! 2 



4tt 2 



/' 



4tt 2 ' 



B, 



(1) 
(2) 
(3) 



which corresponds to the simple pressure-density relation 

1 



P = -(e-4B) 



(4) 



The value of B is not known; however, some constraints 
can be derived. We can specify a lower limit of B due 
to the fact that nucleons do not decay spontaneously to 
two-flavor quark matter. Madsen [25[ shows that this 
lower limit is B 1 / 4 > 145 MeV and gives an expression 
for the energy per baryon E/A as function of B, 



E/A = 829 MeV 



B 1 ! 4 
145 MeV 



(5) 



Since nuclear matter has an energy per baryon of E/A ~ 
930 MeV, according to §5§ bag constants lower than 
B 1 / 4 = 160 MeV correspond to absolutely stable SQM. 

The next step to a more realistic EoS would be to in- 
clude the masses of the quarks. Although the current 
masses of u- and d-quarks are at most 10 MeV and are 
therefore negligible, the mass of the s-quark is of the 
order of 100 MeV. However, in this case an analytic ex- 
pression for P, e and n is no longer possible for finite 
temperatures. Including quark masses as well as QCD 
interactions [26[ leads, for example, at a given B to a the 
energy per baryon which is about 20 MeV higher than 
given by © 27] and thus shifts the range of bag con- 
stants in which SQM is absolute stable. 
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III. COMBUSTION 

We model the conversion from hadronic matter into 
SQM as a combustion, initiated by a seeding of SQM 
which we assume to occur in the center of the star. We 
do not specify the origin of the initial SQM seed (see [|[ 
for various possibilities, or for new ideas). The flame 
front, initially consisting of the boundary surface of some 
central seed, propagates outwards and converts hadronic 
matter into SQM, provided this reaction is exothermic. 
If this is the case, the difference in the energy per baryon 
is released into internal energy and therefore the temper- 
ature increases. The analogous case in chemical combus- 
tion theory is called premixed combustion, where fuel and 
oxidizer are already mixed at low temperatures and the 
flame propagates by conduction of heat [28| . In the case 
of the burning of hadronic matter into SQM the abun- 
dance of s-quarks plays the role of temperature; accord- 
ingly the diffusion of s-quarks leads to the propagation 
of the flame front. The combustion process takes place 
on length scales of the micro-physical reactions, which 
can be estimated as follows: The disintegration of a nu- 
cleon into quarks happens on time scales of the strong 
interaction, ~ 10 _24 s, corresponding to a length scale of 
~ 1CP 13 cm. The conversion of a d-quark into an s-quark 
due to the weak interaction takes place in ~ 10 -8 s. Since 
the weak processes are much slower, they determine the 
time scale of the burning, leading to a width of the reac- 
tion zone, Zbum, not exceeding 10 2 cm, whereas realistic 
calculations yield ibum ~ 10 cm These length scales 
are much smaller than the resolution we can achieve in 
our simulations (Resolved > 10 3 cm) and therefore we can- 
not resolve the reaction zone. Instead, we model the 
conversion front as a discontinuity which separates the 
"unburnt" (hadronic) matter from the "burnt" (strange 
quark) matter and have to take the propagation velocity 
of the conversion front with respect to the fluid flow as an 
input parameter, since this velocity is not determined by 
the hydrodynamic equations but by micro-physical pro- 
cesses on scales of the internal structure of the conversion 
front. 

A combustion can take place either as a supersonic 
detonation driven by a shock wave, or as a subsonic de- 
flagration driven by diffusion processes. Since we can- 
not resolve the internal structure of the flame we have 
to decide before starting our computations whether to 
model the conversion as a deflagration or as a detonation. 
Drago et al. [l2j examine the conversion of hadronic mat- 
ter into quark matter based on the hydrodynamic jump 
conditions. They assume the combustion to start as a 
deflagration and conclude that the process should s tay 
subsonic. Also Niebergal et al. [13| and Horvath [l6[ 
assume the conversion to be subsonic. Based on these 
recent publications we decided to choose a deflagration 
as combustion mode, though we do not exclude the det- 
onation mode and might consider it in future work. 

The relevant input velocity for a deflagration is the 
laminar burning velocity i>i am , which is only very poorly 



known for the burning of hadronic matter into SQM. 
The first attempts to determine it where made by Olinto 
[l4| , who estimates i>i am based on the diffusion of strange 
quarks and the equilibration of the SQM via weak inter- 
actions. The resulting velocities are generally rather low 
but strongly temperature dependent and would lead to a 
wide range of neutron star conversion timescales from 
milliseconds up to several minutes. Recently, Nieber- 
gal et al. [IH conducted one-dimensional hydrodynamic 
simulations of the combustion flame, including neutrino 
emission and strange quark diffusion. They found lami- 
nar burning velocities much higher than in earlier work. 
Because the methods of Niebergal et al. [l3| are more 
sophisticated than in previous publications, we adopt a 
weakly density-dependent laminar burning velocity based 
on a linear fit to their results. This leads to vi am ~ 
10 8 cm/s in the center of the initial neutron star at densi- 
ties of e ~ 10 15 g/cm 3 . Since according to our simulations 
the burning velocity is strongly enhanced by turbulence, 
the importance to know the exact value of vi am is rather 
subordinate (see below). 



A. Turbulent Combustion 

Under certain conditions the laminar propagation of 
the conversion front can be distorted by Ray leigh- Taylor 
(buoyancy) instabilities [see [2i| and references therein] . 
A necessary condition for this is that the gradient of the 
gravitational potential and the gradient of the total en- 
ergy density point in opposite directions ( "inverse density 
stratification" ) . 

In chemical flames, as well as during the thermonu- 
clear burning of carbon and oxygen in the center of a 
white dwarf, the large amount of energy released during 
the burning process leads to a sharp increase in tem- 
perature. In chemical flames a strong increase of pres- 
sure, or a strong decrease in density at constant pres- 
sure, is the natural result and therefore is usually taken 
for granted in qualitative considerations. Similarly, in 
SNe la the degeneracy of the matter is partially lifted, 
therefore the density decreases also in this case, albeit 
not as strongly as in chemical flames. Moreover, in these 
cases, although the chemical abundances change during 
the burning process, the EoS does not change dramati- 
cally. In the case of the burning in white dwarfs at den- 
sities ^ 7-8 x 10 9 g/cm 3 this leads to an inverse den- 
sity stratification, instabilities and turbulence [29| . How- 
ever, because of the strongly degenerate state of matter 
in neutron stars and the fundamentally different EoS be- 
fore and after the conversion process it cannot be taken 
for granted that the neutron star matter behaves in the 
same way as described above. The state of the fluid be- 
hind the conversion front is determined by the change 
of the EoS and the hydrodynamic jump conditions [see 
e.g. LL2( which result from the conservation of the baryon 
flux density and the energy-momentum tensor across the 
flame surface and has to be computed in hydrodynamic 
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simulations. To explore if in the vicinity of the propa- 
gation front the density of the SQM is lower than the 
density of the hadronic phase for our choice of EoS is 
therefore one aim of this work. 

The Rayleigh- Taylor instability can only grow and lead 
to turbulent motion if the perturbations of the front ex- 
ceed some minimal length scale, A m i n , which depends on 
the burning velocity, the gravitational acceleration g, and 
the density contrast between the total energy density of 
the hadronic phase eh and the total energy density of the 
quark phase e q [2{| , 

A min = 2nvf am ( g 6h q ) . (6) 
\ £h + e q J 

We calculate A m ; n for different setups in Section IY Al 

In the established heuristic turbulence model [10, HH 
instabilities like the Rayleigh- Taylor instability (and sec- 
ondary shear instabilities) lead to turbulent eddies on 
large scales, which decay successively into ever smaller 
eddies until, at the Kolmogorov length scale Ik, viscosity 
effects dissipate the smallest eddies into thermal energy. 
In this turbulent cascade kinetic energy is transported 
from the largest to the smallest scales and is finally dissi- 
pated. This picture assumes that magnetic fields do not 
significantly affect the dynamics. For the velocity fluc- 
tuation v(l) on a given scale /, which can be interpreted 
as the turnover velocity of an eddy of size ', this model 
yields the Kolmogorov scaling [32| . 

v(l) = v(L) ^j-j , (7) 

where L is the integral scale, the size of the largest eddies. 

The Reynolds number Re on different scales is there- 
fore 

( z\ 4/3 

Re(l) = Re(L) ( -J , (8) 

since Re(l) oc v(l)l. Horvath and Benvenuto @ esti- 
mate the Reynolds number of flows in both neutron and 
strange stars to be Re(L) ~ 10 10 . At the Kolmogorov 
scale Re(l]i) ~ 1 holds, so we get 

and hence Ik « 'bum- 

The scale on which the eddy turnover velocity is equal 
to the laminar burning velocity is defined as the Gibson 
scale l G [e.g.EH, 

v(l G ) = wi am . (10) 

Turbulence cannot distort the flame front on scales 
smaller than Iq since according to (O on these scales 
the eddy turnover velocity is smaller than the laminar 
burning velocity, whereas on scales larger than Iq, the 
turbulent eddies alter the the shape of the flame front. 



If we assume the above scaling law for a ris- 
ing Rayleigh- Taylor unstable bubble of typical size 
L w 10 5 cm and typical macroscopic velocity variations 
v(L) fx 10 9 cm/s, we find Iq = 10 2 cm. The combus- 
tion theory was developed for chemical flames and only 
adapted to SNe la [33l . H3| , whereas no detailed studies 
were conducted for the case treated in this work. How- 
ever, the Kolmogorov scaling was found to fit quite well 
in the case of SNe la [H, [36| , so based on these results 
and in the absence of exact calculations we assume that 
this is the case for our problem as well and obtain 

'burn < Ig < 'resolved- (H) 

This leads to two important consequences: 'burn < 'g 
means that the turbulent eddies cannot disturb the flame 
front. Thus, it can still be described as a well-defined 
discontinuity. The burning is said to take place in the 
flamelet regime [28j]: Although the internal flame struc- 
ture is not disturbed, the total burning rate is enhanced 
as turbulence alters the geometry and thus enlarges the 
surface of the front. Since Ig < 'resolved the surface of the 
flame front is also enhanced on unresolved scales, leading 
to an increase in the effective front propagation velocity 
on these scales. This effective velocity is described by 
the turbulent burning velocity t>t U rb, which is defined as 
the mean propagation velocity of the flame front at the 
marginally resolved scale. 

For strong turbulence, the turbulent burning velocity 
becomes independent of the laminar burning velocity, as 
is the case during the thermonuclear burning of a white 
dwarf. In this work we aim to explore if the same is true 
in the conversion process of a neutron star. 

B. Conditions for Exothermic Combustion 

Since we describe the conversion of hadronic matter 
into SQM as a combustion, and a combustion has to be, 
by definition, exothermic [37J , we can specify the follow- 
ing necessary condition for the conversion to take place: 
The total energy density of the quark phase e q in a ther- 
modynamic state (P, X) has to be lower than the energy 
for the hadronic matter eh in the same state [37| . 

e h (P,X) > e q (P,X), (12) 

where P is the pressure, X is the generalized volume, 
X = (e + P)/n B , and n B is the baryon density. In the 
case of our analytic EoS for SQM (Q}, this can be rewrit- 
ten as a simple condition for the energy density of the 
hadronic phase [t| HH : 

e h {P) > 'iP + AB. (13) 

From this relation it becomes clear that for each given 
total energy density eh and temperature Th the corre- 
sponding pressure of the hadronic phase P and the value 
of B determine whether the conversion can proceed in 
form of a combustion wave. Thus, after choosing the 
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FIG. 1. Maximum bag constant B allowing an exothermic 
combustion as function of the total energy density e for two 
different hadronic EoS (Lattimer-Swesty with K = 180 MeV 
and Shen). The horizontal line indicates the theoretical lower 
limit of B. For each EoS, two cases are plotted: in the first 
case temperature Th and proton fraction Y p are kept constant, 
in the second case the matter is in /3-equilibrium at zero tem- 
perature. 



EoS and assuming a fixed Th we can calculate for each 
eh a critical bag constant, -B C rit(e/i), which is the largest 
possible bag constant for an exothermic combustion. The 
results of these calculations using both the LS EoS with 
K = 180 MeV and the Shen EoS are shown in Figure [1] 
Here the results arc plotted for two different cases: In the 
first case we assume a constant temperature of the un- 
burnt hadronic matter of Th ~ 100 keV and a constant 
proton fraction of Y p = 0.2. We adopt these assump- 
tions for our numerical simulations presented in Section 
IIIII In the second case we assume /3-equilibrium and zero 
temperature. As visible in Figure [TJ the differences be- 
tween the two cases are rather small and thus negligible 
for the qualitative treatment in this work. Also apparent 
from this figure is that for bag constants larger than the 
theoretical lower limit, B 1 / 4 > 145 MeV, and tempera- 
tures found in the interior of cold neutron stars, hadronic 
matter described by the Shen EoS cannot be burned into 
SQM in an exothermic combustion, regardless of the den- 
sity. In contrast, matter described by the LS EoS can be 
converted into SQM in an exothermic way at densities 
occurring in the center of neutron stars. The difference 
between the two hadronic EoS can be explained as fol- 
lows. The Shen EoS is rather stiff, much stiffer than the 
LS EoS, that is at the same density the pressure is much 
higher. According to ([TT?)) this leads to a higher energy 
threshold for a given density. Based on this results we 
have to refrain from using the Shen EoS in our hydrody- 
namic simulations. 

The LS EoS can be used with different incompress- 
ibility moduli K, we consider K — 180 MeV and K = 



FIG. 2. Maximum bag constant B allowing an exothermic 
combustion as function of the total energy density e for the 
LS EoS and two different incompressibility moduli K. The 
horizontal line indicates the theoretical lower limit of B. 
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FIG. 3. Minimum total energy density e for an exothermic 
combustion as function of the bag constant B and for different 
temperatures Th, using the Lattimer-Swesty EoS with K — 
180 MeV. 



220 MeV. We compare these two possibilities in Figure 
[2] For low bag constants (B 1 / 4 - 145 MeV) the higher 
stiffness of the EoS with higher K affects the lower den- 
sity limit only slightly, but for B 1 / 4 > 152 MeV the range 
in which exothermic combustion is possible becomes very 
narrow. Since our goal is to conduct simulations with 
higher bag constants to be able to compare the results 
for a wide range in the amount of released energy, we use 
in our simulations only the LS EoS with K = 180 MeV. 
In Figure [3] we concentrate on this case. Here we plot 
the minimum total energy density of the hadronic phase, 
e m in(-B), as a function of B and for different fixed tem- 
peratures. Since below this density threshold no combus- 
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tion is possible, it plays an important role in our simu- 
lations. The continuous line in Figure [3] shows the case 
with Th = 100 keV, the temperature we adopt for the cold 
neutron star in our simulations. In addition we explore 
the effects of several higher temperatures. For temper- 
atures up to Th = 1 MeV only slight differences would 
be visible due to the strong degeneracy of the matter. 
In proto-neutron stars considerably higher temperatures 
occur, therefore also results for = 10 , 20 and 30 MeV 
are shown in the figure. These temperatures have a no- 
ticeable effect on the density threshold, as visible in Fig- 
ure [3J In general, higher temperatures move the density 
threshold to higher densities. 



IV. NUMERICAL METHOD 

Taking advantage of the methodical similarities of com- 
bustion in white dwarfs and in neutron stars, we use an 
existing code, which is well tested and frequently applied 
for various SN la related simulations [e.g. HH, [401 an d 
adapt it for the subject of this work. The main features 
of this code will now be described briefly. 

The reactive Euler equations are solved using an ex- 
plicit piecewise parabolic method (PPM) (4lj . a higher- 
order Godunov scheme. Specifically, our code is based 
on the Prometheus implementation [42[ of PPM. 

To track the flame front we use the level-set method 
which was introduced by Osher and Sethian I4p and im- 
plemented in the code by Reinecke et al. [4J|. In this 
scheme, a signed distance function G, which is positive 
in the burnt material and negative in the unburnt mate- 
rial, is assigned to each point in the computational do- 
main. The zero level-set of G thus separates the burnt 
from the unburnt matter and marks the location of the 
flame front. The level-set is propagated with the burning 
velocity perpendicular to the flame surface and advected 
as a passive scalar without fundamental modifications of 
the hydrodynamics solver. It is now possible to calculate 
the burnt and unburnt volume fractions in each cell. 

To ensure that the regions of highest interest are op- 
timally resolved for a given fixed number of grid cells 
and no computational resources are wasted on regions 
of subordinate importance, the computational domain is 
separated into two grids [45|, H(| , an outer coarser grid, 
where the cell size increases outwards, and an uniformly 
spaced moving inner grid tracking the conversion front 
and expanding with it into the outer grid. This way we 
achieve an initial resolution in the center of the star of 
~ 2.6 x 10 3 cm x (grid cells per dimension/128) -1 , if our 
grid covers one octant of the star. 

As described in Section MI A| we cannot resolve the 
turbulent motion down to the Gibson scale. Therefore, 
we perform "large eddy simulations": Only the largest 
scales of the system are resolved, while the turbulent 
motion on smaller scales is modeled by a subgrid scale 
(SGS) model. The SGS model determines the turbulent 
energy, from which the turbulent burning velocity can 



be inferred. Schmidt et al. [13, |48j] introduced a sophisti- 
cated localised SGS turbulence model and implemented 
it into the code. This model determines the SGS turbu- 
lence velocity <?sgs • The turbulent burning velocity wturb 
is then obtained by setting [48| 



wt - b=uiam v 1+ K?f) ' (14) 

with the laminar burning velocity wi am as a lower limit. 

The code as described was written to set up a white 
dwarf in a Newtonian gravitational potential and to 
model the thermonuclear burning of carbon and oxygen. 
Thus, it had to be adapted to the subject of this work. 
In contrast to the case of white dwarfs (compactness 
(GM/_Rc 2 )wd ~ 0.001), in neutron stars (compactness 
(GM/_Rc 2 )ns ~ 0.1) general relativistic effects cannot be 
neglected. Computations in full general relativity are, 
however, beyond our scope. Given the overall uncertain- 
ties, particularly in the EoS, we consider the error in- 
troduced by the use of Newtonian dynamics to be not 
critical, however a comparison of our results with gen- 
eral relativistic simulations would be interesting. But 
a modification of the gravitational potential cannot be 
avoided, otherwise the results would be completely be- 
side the point. For example for a given mass of the neu- 
tron star the central density would be much lower and 
thus exothermic combustion would not be possible at all. 
Therefore an effective relativistic gravitational potential 
[49j based on the Tolman-Oppenheimer-Volkov (TOV) 
equations was implemented. 

Further adaptions of the code to the new setup include 
the EoS for hadronic and quark matter as described in 
Section UH and, as a replacement for the "burning rou- 
tine" in the original code, a routine which takes care of 
the conversion of the hadronic matter into SQM. This 
takes place by switching to the quark EoS and releas- 
ing the difference in the energy per baryon into internal 
energy, while conserving the total energy. 

We set up one octant of the neutron star on a three- 
dimensional Cartesian grid with 128 or 192 grid cells in 
each dimension and applied reflecting boundary condi- 
tions at all borders of the computational domain. Burn- 
ing is initialized in the following way: At the center of 
the star we construct a small sphere with a radius of 
fsccd = 10 5 cm on which a sinusoidal perturbation with 
an amplitude of 2 x 10 4 cm is superimposed. The initial 
seed is shown in the close-up of Figure [5] (a). When start- 
ing the simulation, the matter inside this small volume is 
converted instantly and constitutes the initial SQM seed. 

Since both the size and the form of the initial seed 
are not known, we choose this configuration for numer- 
ical reasons: The size of the perturbations is similar to 
the minimum length scale for turbulent burning A m i n (cf. 
Sections llII Al and IV Al) . therefore the front is expected to 
develop Rayleigh- Taylor instabilities soon after the start 
of the simulations. Smaller initial perturbations would 
need some time to grow before Rayleigh- Taylor instabili- 
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Model 


Resolution 


B 1/4 /MeV 


M un burnt / Mq 


B147.128 


128 3 


147 


0.48 


B150.128 


128 3 


150 


0.66 


B150.192 


192 3 


150 


0.67 


B152.128 


128 3 


152 


0.77 


B155.128 


128 3 


155 


0.99 



TABLE I. Overview of the different models. M U nburnt is the 
gravitational mass of the remaining hadronic outer layer at 
t = 3.0 ms, when the combustion can be considered as com- 
plete in all cases. 



ties become possible. But since in the end the core is con- 
verted completely, the results should change only slightly, 
whereas the computational costs would be considerably 
higher. As described in Section [ill! we assume the com- 
bustion to be a deflagration and ignite the burning ac- 
cordingly. Although we do not expect different initial 
configurations to alter our results considerably, possible 
effects of different initial geometries and different ways of 
ignition will be explored in future work. 




t [ms] 



FIG. 4. Gravitational mass of unburnt (hadronic) material 
in the three-dimensional simulations for different bag con- 
stants B, as a function of time (models B155_128, B152_128, 
B150_128 and B147_128). 



V. SIMULATIONS 

We conduct several runs with varying bag constant B. 
Since only some constraints on B are known, we can use 
it as a parameter to change the EoS for SQM and are 
thus able to control the amount of released energy from 
very high to rather low values, cf. ([5]). We vary B in 
a subset of the theoretically admissible range between a 
lower limit of = 147 MeV and an upper limit of 

^Wgh = 155 MeV. At even higher B, the combustion 
would be restricted to the very innermost region of the 
neutron star or would not be possible at all, cf. Fig- 
ure CD and CE!]). We use B 1 / 4 > 155 MeV only to test 
if instabilities grow at the beginning of the burning; re- 
sults are presented in Section IV Al In alternative units 
our chosen limits are roughly B\ ow — 60MeV/fm 3 and 
£?high = 80MeV/fm 3 - values also used as a lower and 
upper limits in the literature [e.g. [27| . 

We start our computations with an non-rotating, cold, 
isothermal "standard neutron star" in hydrostatic equi- 
librium, having an initial central total energy density of 
e c = 1.0 x 10 15 g/cm 3 , a gravitational mass of M — 
1.4M , a radius of R = 11km, a proton fraction of 
Y p = 0.2, and a temperature of T = 100 keV. We con- 
ducted four runs with a resolution of 128 grid cells per 
dimension and bag constants of B 1 / 4 = 147 , 150 , 152 
and 155 MeV, respectively. Table U shows an overview 
of the models presented here. In Figure |4] the tempo- 
ral evolution of the conversion for different B is shown, 
represented by the gravitational mass of the remaining 
unburnt hadronic material. 

In addition, we conducted one run with a higher resolu- 
tion, 192 grid cells per dimension, using an intermediate 




0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 

t [ms] 



FIG. 5. Resolution study: Two models with B 1/4 = 150 MeV 
which differ only in resolution (B150_192 and B150_128). 



bag constant of B 1 / 4 = 150 MeV (model B150.192). To 
study the effects of different resolutions, we compare in 
Figure [5] the two models B150.128 and B150_192, which 
differ only in the resolution (128 3 and 192 3 , respectively). 
Apparently there are only slight differences between the 
two models. In particular the slopes in the phase of rapid 
burning, which are determined by the conversion rate, 
which in turn depends on the turbulent burning velocity, 
agree very well. The different resolutions only become 
noticeable in the representation of the exact position of 
the density threshold for exothermic combustion - hence 
the slight discrepancy in the amount of unburnt matter at 
later times. Therefore we consider our simulations con- 
verged in the sense that the effects caused by resolution 
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B 1/4 /MeV 


145 


147 


150 


152 


155 


157 


At 


0.11 


0.091 


0.067 


0.051 


0.027 


0.010 


A m in/10 4 CHI 


3.6 


4.4 


6.2 


8.5 


16 


45 



TABLE II. Atwood number At and minimal length scale for 
turbulent burning A m i n for different bag constants B at time 
t = 0.1ms, determined in three-dimensional simulations. 




10 11 12 13 14 15 16 

distance from center [10 4 cml 



17 



FIG. 6. Comparison of the minimal length scale for turbulent 
burning A m i n in the early phase of the conversion process for 
different bag constants B and points in time, determined in 
three-dimensional simulations. The number on each line in- 
dicates B 1 / 4 in MeV. For each B the first and second point 
correspond to time t\ = 0.1ms and t^ = 0.2 ms, respectively. 
On the abscissa the average position of the conversion front 
at ti and ti is shown. 



are smaller than uncertainties caused by other sources. 
Thus, we regard a resolution of 128 cells per dimension 
to be sufficient for our quantitative analysis. 

After addressing the question of whether burning is 
turbulent, the results of the simulation with the highest 
resolution, model B15CL192, are discussed in some detail 
below. In the subsequent sections we will briefly discuss 
differences in the two extreme cases (models B147_128 
and B155_128). 



10 1 



i<r 



10° 



10 





maximum burning velocity 

— average burning velocity 

average laminar burning velocity 


) J 







0.0 0.5 1.0 1.5 



2.0 2.5 

t [ms] 



3.0 3.5 4.0 4.5 



FIG. 7. Burning velocity: Comparison at each timestep of 
maximum burning velocity, average burning velocity and the 
underlying average laminar burning velocity. The averages 
where done over all cells in which burning occurs. Data 
from the high resolution run with B 1/A = 150 MeV (model 
B150_192). 



termincd at t = 0.1ms. As visible in Table ITTI and Fig- 
ure [6) A m in depends strongly on B, and becomes very 
large for high B. For the highest examined bag constant, 
B 1 / 4 = 157 MeV, A m j n is comparable to the size of the 
system and no growth of Rayleigh- Taylor instabilities is 
expected. Bag constants starting at B 4 ! 4 — 152 MeV 
down to the lowest B, lead to smaller A m i n which are 
comparable to or smaller than the size of the initial per- 
turbations - thus instabilities can grow. In the simulation 
with B 1 ! 4 — 155 MeV this is not the case at t = 0.1ms 
but already at some slightly later time, since A m i n de- 
creases with time as the gravitational acceleration be- 
comes stronger, cf. ([6|) and Figure [6] Our simulations 
confirm this: in all runs except for B 1 / 4 > 157 MeV we 
see Rayleigh- Taylor instabilities form. Thus the burning 
of a neutron star into a quark star becomes turbulent in 
most cases, given our choice of EoS. 



Intermediate Case: B 1/4 = 150 MeV 



A. Onset of Turbulence 

We calculate the minimum length scale for turbulent 
burning, A m j n , according to ([B]), see Section IlII Al To 
ensure comparability, we use the same three-dimensional 
setup for all B, as described above, and the same resolu- 
tion of 128 grid cells in each dimension. 

In Figure [5] we compare A m j n at the beginning of the 
conversion process for different bag constants B and 
points in time. The density contrast is quantified by 
the Atwood number At — (e/j — e q )/(eh + e q ). Table [TT1 
lists At and A m ; n for different B. The values were dc- 



In this section we present a detailed discussion of the 
results of the simulation with a resolution of 192 grid 
cells per dimension and an intermediate bag constant, 
B 1 / 4 = 150 MeV (model B150.192). According to © 
the energy per baryon in this case is E/A — 858 MeV, 
corresponding to a difference of ~ 70 MeV per baryon 
with respect to the energy of nuclear matter. 

In Figure [5] (a) the initial configuration including the 
SQM seed in the center can be seen. The shape of the 
seed as described in Section llVI is shown additionally in 
the close-up in this figure. After ignition the conversion 
front propagates into the hadronic matter, at first in a 
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0.2 0.4 0.6 0.8 1.0 1.2 1.4 



(a) t = 



(b) t = 0.7 ms 





0.2 0.4 



1.0 1.2 1.4 



(c) t = 1.2 ms 



(d) t = 4.0 ms 



FIG. 8. (color online) Model B150_192: Conversion front (red) and surface of the neutron star (yellow) at different times t. In 
(a) a close-up of the central region is added. Spatial units 10 6 cm. 



laminar way until initial perturbations of the conversion 
front become unstable due to Rayleigh- Taylor instabili- 
ties. Until turbulence has fully developed, the conversion 
process stays in a short phase of nearly laminar burning 
while the instabilities grow, see Figure U which shows 
the amount of unburnt (hadronic) matter as a function 
of time, and Figure[7j where we compare the average lam- 
inar burning velocity, the average burning velocity and 
the maximum burning velocity at each timestep. The av- 
eraging was done over all cells in which burning occurs. 

As the instabilities grow, typical mushroom-shaped 



structures, rising plumes of SQM, are forming and 
hadronic matter is falling down in between. These struc- 
tures can be seen in Figure [5] (b) , where the conversion 
surface is shown at t — 0.7 ms. Starting at t ~ 0.5 ms 
strong turbulence develops and rapid burning takes place 
until t r~ 1.5 ms, as visible in Figure [4] The structure of 
the conversion front near the end of this phase of rapid 
burning can be seen in Figured] (c) . The plumes grow un- 
til the conversion front reaches densities where the condi- 
tion for exothermic combustion (1121) is no longer fulfilled. 
They continue to grow laterally, until they eventually 
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merge, leaving bubbles of hadronic matter in between. 
Turbulence then weakens and the flame slows down. The 
remaining pockets filled with hadronic material shrink 
until they eventually vanish completely. Now all matter 
inside the volume confined by the above mentioned den- 
sity threshold is burned and the star consists of an inner 
sphere of SQM containing about half of the mass and an 
outer layer of unburnt hadronic matter (cf. Figure [5] (d)). 
This outer layer has a mass of about 0.67 M Q and densi- 
ties lower than the threshold (fT2"j) but mostly still super- 
nuclear (applying LS EoS and B 1/4 = 150 McV, the den- 
sity threshold in cold matter is at about 1.8 p nuc i car , see 
Figure [3]). 

Turbulent motions leads to burning velocities consid- 
erably higher than the laminar burning velocities, the 
maximum amplification factor is about 50 and on aver- 
age between 2 and 20 (see Figure [7]). This figure also 
clearly shows that the turbulent burning velocity and 
thus the strength of the turbulence increases rapidly un- 
til it reaches a maximum at t ~ 1.0 ms. At that point a 
steady but slower decrease starts. The maximum Mach 
numbers reached were about 0.2, i.e. the combustion was 
clearly subsonic. As we do not include any kind of cool- 
ing, the large amount of energy released in the burning 
process is turned into thermal energy and the inner SQM 
region is heated to temperatures of about 50 MeV in the 
center of the star. 

We stopped this simulation at t = 4.0 ms. By then 
the conversion rate has dropped to a very low value and 
seems to approach zero asymptotically. Since at that 
time the system is approximately in hydrostatic equi- 
librium (the dynamical time scale of a neutron star is 
Tdyn ~ 5 x 10 -2 ms) we do not expect any further con- 
version of a significant amount of mass. Therefore the 
structure of the remnant should not change if the simu- 
lation would have been carried on for longer timescales 
- at least in our model without cooling processes and in 
the approximation of a hydrodynamic combustion. 

C. Lower Limit: B^* = 147 MeV 

Now we briefly discuss the simulation with 128 grid 
cells per dimension and with our lower limit for the bag 
constant, B^ = 147 MeV (model B147_128). This cor- 
responds to the largest difference in energy per baryon 
compared to nuclear matter, E/A = 90 MeV. 

Qualitatively, the conversion process evolves in the 
same way as in the case described above (model 
B150_192), but there are some quantitative differences: 
The energy release is higher than in the intermediate 
case, therefore the burning leads to a stronger inverse 
density stratification, resulting in a faster growth of in- 
stabilities and stronger turbulence. The rising plumes of 
SQM are observable in Figure [9] (a) as typical "mush- 
rooms" , like in the previous case. Comparing Figure [5] 
(a) and Figure [5] (b), both showing the conversion surface 
at t — 0.7 ms, clarifies that the conversion process takes 



place considerably faster for the lower B. Figure [4] shows 
that after a short phase of slow burning, rapid burn- 
ing occurs from t ~ 0.4 ms until t ~ 1.5 ms. Then the 
burning slows down and the conversion rate approaches 
zero. At t = 5 ms, the remnant has an inner SQM core 
with a radius of ~ 9km, cf. Figured] (b), surrounded 
by an hadronic outer layer with a mass of 0.48 Mq, the 
least massive outer layer in all our simulations. Central 
temperatures of the core reach 53 MeV, somewhat higher 
than in the previous case due to the higher energy release. 



D. Upper Limit: B^ h — 155 MeV 

Finally we present the simulation with 128 grid cells 
per dimension and our highest bag constant, -B^gh = 
155 MeV (model B155.128). Here the difference in en- 
ergy per baryon, E/A ~ 40 MeV, is considerably lower 
than in the cases B147.128 and B150.128. Figures [TO] 
(a) and [T0l (b) show the conversion front at t — 0.7 ms 
and at the point when we stopped our simulation, at 
t = 4.6 ms. From the figures the similar evolution com- 
pared to the above described cases with lower B are vis- 
ible. The lower E/A and the higher density threshold 
for exothermic burning (cf. Figure [3|) lead to a slower 
and less violent burning, which ceases at higher densities 
compared to the models previously shown. Consequently, 
at the end of the simulation the resulting strange matter 
core is smaller and is surrounded by an hadronic outer 
layer of 0.98 M Q . Temperatures of around 45 MeV are 
reached in the center. Figure [4] shows that the conver- 
sion rate, represented by the slope of the curves, is lower 
than in the other cases and the combustion takes longer 
although less material is burnt. 



VI. CONCLUSIONS 

We presented three-dimensional hydrodynamic simu- 
lations of the conversion of a neutron star into a quark 
star assuming different bag constants B for describing 
SQM. In all cases we observe growing Ray leigh- Taylor 
instabilities of the conversion front. The resulting tur- 
bulent motion enhances the conversion velocity strongly, 
leading to conversion timescales of Tbum ~ 2 ms for all B. 
However, recent suggestions [H, HH that the turbulence 
enhances the burning speed to sonic or even supersonic 
velocities could not be confirmed, which came as no sur- 
prise since in the analogous case of SN la such a transition 
is not possible either as long as burning proceeds in the 
flamelet regime |33j |. 

In all cases we observe at the end of our simulations 
a spherical SQM interior surrounded by an outer layer 
of hadronic matter. This outer layer exists because in 
our hydrodynamic approximation the combustion stops 
when the conversion front reaches conditions under which 
exothermic burning is no longer possible. Since this 
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(a) t = 0.7 ms 

FIG. 10. (color online) Model B155.128: Conversion front (rei 
Spatial units 10 6 cm. 

condition depends on density and is fulfilled for suffi- 
ciently high densities only, it can roughly be described as 
a density threshold which forms a boundary that sepa- 
rates the high density (burnt) strange quark matter and 
the low density (unburnt) hadronic matter. In our ap- 
proximation we can make no statement on whether the 
conversion process proceeds further beyond this bound- 
ary by processes which cannot be described as a com- 




(b) t = 4.6 ms 

and surface of the neutron star (yellow) at different times t. 



bustion. Possibly free neutrons diffuse into the quark 
matter and are converted subsequently (T3 |. a process 
that probably is exothermic, as Lugones et al. @ already 
pointed out. Free neutrons are abundant in hadronic 
matter at densities higher than the neutron drip density, 
e dri P ~ 4x 10 11 g/cm 3 . However we expect these ad- 
ditional processes to happen on much longer timescales 
than the combustion described in this work. 
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The obvious consequence of an at least temporary ex- 
istence of an outer layer of unburnt hadronic matter is 
that the resulting quark star could support a rather thick 
crust, unlike bare strange stars, which can presumably 
support only a tiny crust. This would allow e.g. for pul- 
sar glitches, if the timescale of the conversion after the 
combustion has ceased is large enough. Some authors 
suggested that the conversion of a neutron star into a 
strange star may eject neutron-rich material from the 
surface, and that in this ejecta the nucleosynthesis of 
heavy neutron-rich nuclei via the r-process may occur 
[50j]. However, our results suggest that ejection of matter 
from the star is rather unlikely since the violent burning 
ceases before reaching the surface. Any subsequent con- 
tinuation of the conversion by processes not describablc 
by a combustion is expected to be much slower, and to 
take place in a much less violent way. But given our ig- 
norance about these processes more detailed work on this 
subject may lead to differing conclusions. 

The existence of the hadronic outer layers, or the possi- 
bility of exothermic combustion even in the center of neu- 
tron stars, depends (like many other properties) strongly 
on the EoS used for the hadronic as well as for the quark 
phase. Hence any firm prediction needs a more realistic 
treatment. Furthermore, the maximum mass configura- 
tion of non-rotating stars of both the LS EoS and our 
bag model EoS have M max < 2 M Q and therefore conflict 
with observations [22j. As mentioned before, we never- 
theless use those EoS in this work because we consider 
them as sufficient for our first attempts. In future work 
we want to improve on this and plan to use more realis- 
tic EoS. Regarding the quark phase, finite strange quark 
masses and QCD-interactions can be included into the 
bag model. SQM bag model EoS which contain these cor- 
rections can support a 2 Mq neutron star, as was shown 



by Weissenborn et al. 51] . Recently also the choice of 
micro-physical finite temperature EoS for nuclear matter 
has become larger [e.g. [52], [53[, so we can consider ad- 
ditional hadronic EoS apart from the LS and Shen EoS 
which we used in this work. Another possibility is to con- 
sider the use of modern zero-temperature micro-physical 
EoS together with an ideal gas component to account for 
temperature effects, whose reliability has been tested in 
54| . Further improvement would be achieved by adding 
neutrino cooling, which could be relevant since rather 
high temperatures are reached in the quark core. Until 
now we use an initial model resembling an old isolated 
neutron star, the same calculations could be done with 
a young (proto)neutron star and in connection with a 
core collapse supernova. Furthermore our computations 
should be extended to make statements about observ- 
able quantities. Therefore we plan to calculate the grav- 
itational wave signal of the conversion of a neutron star 
into a quark star. 
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